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ABSTRACT 



In the present study, a three-dimensional differential field model based 
on general orthogonal coordinate systems is developed. The model, which is 
essentially designed for closed spaces, also includes the physical effects of 
turbulence, strong buoyancy, full coinpressibi 1 ity, pressure rise due to fire 
loading, surface-surface and surface-flame radiation exchange, and heat losses 
through the wall. It is based on a control -vol ume staggered-celi finite- 
difference approach with primitive variables. Results of numerical calcula- 
tions based on the field model are compared with test data for a methanol fire 
in the NRL FIRE I test facility which is in the form of a closed pressure 
vessel. Reasonable comparisons of the resulting pressure and temperatures at 
several locations have been obtained. Also shown are the detailed velocity 
and temperature fields inside the vessel at different time instants after the 
commencement of the fire. 
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INTRODUCTION 



It is now generally recognized that fire hazards in rooms, compartments, 
corridors, passageways, and other confined spaces can only be minimized and 
mitigated by a multi-prong approach involving different strategies. Among them 
is fire modeling which deals with the development of mathematical models which 
describe the physical and chem-ical processes in the spread of fire in spaces as 
a function of the ignition source, space geometry, and material content. Once 
validated by experiments in small-scale laboratory tests and/or full-scale fire 
tests, these models can be used as computer-based simulation models to determine 
the effects of the significant parameters on the fire-spread phenomena. The 
results can then be utilized for developing fire mitigation measures as well as 
for providing a rational basis for post-fire investigations. One chief advan- 
tage of fire modeling is the significantly reduced need for full-scale fire 
tests which are extremely expensive and time consuming. 

Fire models can be categorized as either zone models or field models. Zone 
models deal with dividing the fire-affected environment or space into distinct 
regions or zones which can be separately analyzed either empirically or theore- 
tically in terms of input and output information based on mass and energy balan- 
ces. Examples of such zones include the fire envelope, fire plume, hot gas 
layer, ceiling jets, wall jet, flow through openings, etc. These input and 
output quantities from individual zones are then assembled to describe the 
overall gas dynamic phenomenon of the entire fire environment. In most cases, 
there results a set of ordinary differential equations in time which are then 
solved numerically by standard integration algorithms. While the zone models 
are generally computationally efficient, shortcomings do exist in that the 
models for some of the zones are not adequately known and quantified, and some- 
times there is a question as to exact range of validity of such zone models. 
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Ef forts are now being continued to refine the existing zone models and also to 
compare with full-scale test data for validation purposes (for example see ’ 

Rockett et al., 1987). 

Field models which are based on the numerical solutions to the governing 

I 

differential field equations of the conservation of momentum, mass, energy and 
spices, on the other hand, have not been widely used in the past for fire 
modeling purposes even though these models are inherently more rational and 
capable of revealing both the important large-and small-scale phenomena in a 
given fire-spread problem. The primary reasons are the excessive use of com- 
puting resources and degrees of uncertainty in modeling such phenomena as tur- 
bulence, combustion and thermal radiation. However, particularly because of the 
increasing availability of high-performance computers including supercomputers, 
field models for fire modeling have found increasing use in more recent times. 
For instance, such models have been developed at the National Bureau of 
Standards (Baum and Rehm, 1984), Imperial College of London (Bagnaro, Laouisset 
and Lockwood, 1983), Borehamwood Fire Research Station (Markatos et al . , 1982), 
and the University of Notre Dame (Yang et al., 1984, Yang and Lloyd, 1985), Kou 
et al., 1986, and Nies, 1986). In view of the different physical and chemical 
models utilized in each of these field models and the fact that differences also 
exist in the details of the numerical algorithms involved, there is no way to 
ascertain the relative accuracies of these models. Comparisons of the numerical 
results with full-burn tests are scarce, but do exist and are fairly reasonable 
(Yang et al., 1984; Yang and Lloyd, 1985; Bagnaro et al., 1983, Markatos and 
Pericleous, 1983). An ultimate assessment of the various field models cannot be 
made until more validation studies based on full-scale experimental data become 
available. Additional discussions on fire modeling have been recently given by 



Stroup (1987). 
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Despite the above-mentioned uncertainties of the field models, many of the 
important physical effects can now be accommodated. Those include turbulence 
(Yang and Liu, 1978; Bagnaro et al., 1983; Marketos and Pericleous, 1983), ther- 
mal radiation (Lloyd et al . , 1979; Bagnaro et al., 1983; Markatos and 
Pericleous, 1983); Nies, 1986), combustion (Bagnaro et al., 1983; Markatos and 
Pericleous, 1983), three-dimensional vented spaces (Bagnaro et al., 1983; 

Markatos and Pericleous, 1983; Yang and Lloyd, 1985; Nies, 1986), closed-space in- 
duced pressure rises (Nies, 1986), and wall losses (Nies, 1986). However, all the 
above-mentioned field models deal with box-like rectangular spaces or enclo- 
sures, and yet, many full -burn facilities have other geometries for which these 
field models are very difficult to apply. A good example is the facility known 
as FIRE I located at the Naval Research Laboratory (NRL), which is in the form 
of a pressure vessel including a central cylindrical body plus two spherical end 
caps (Alexander et al., 1982). The purpose of this report is to present another 
field model which is based on a generalized orthogonal coordinate system which 
still retains the same capabilities in terms of physical models for turbulence, 
full compressibility, strong buoyancy, surface-to-surface and f lame-to-surfacs 
■•adiation, wall losses, and effects of heat addition in a closed space. The rec- 
tangular coordinate-system case as given by Nies (1986) becomes a special case 
of this new field model, and it can be mentioned that the model of Nies (1986) 
is basically an extension of the original two-dimensional UNDSAFE-II code devel- 
oped at the University of Notre Dame (Yang and Liu, 1973) for vented rec- 
tangular enclosures. Furthermore, this new field model will be demonstrated ny 
presenting a set of numerical result.s which simulate some experimental data 
obtained in the NRL FIRE I facility. Finally, several additional capabilities 
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of this new field model designed to better simulate real fire scenarios will 
also be delineated. 

MATHEMATICAL FORMULATION OF THE FIELD MODEL 

As pointed out previously, all field models for the fire-spread phenomena 
are based on first-principle conservation differential equations for mass, 
momentum, energy and species. The approach utilized here is similar to that 
given by Yang et al. ( 1987a, b), and basically, the governing equations in the 
Cartesian coordinates are first transformed into those of a generalized orthogo- 
nal coordinate system, which are then f inite-differenced by the control -volume 
staggered-cel 1 primitive-variable approach. The flux terms are handled in such 
a way that the final finite-difference algebraic equations resemble closely 
those of the Cartesian coordinates. As a result, the computer algorithm for the 
generalized orthogonal coordinates is no more complex then that for the 
Cartesian coordinates. In fact, it has recently been shown that the present 
approach represents a viable alternative to the popular approach of utilizing 
body-fitted coordinates (Yang et al., 1988). Under the conditions of turbulent 
compressible flow with buoyancy, the governing equations for the mean quantities 
in the Cartesian coordinates may be written as follows in the standard tensor 
forms: , 

Pt + ( PUi) = 0 (1) 

(pui)t + (PUiUj) J - PGi + (Pij),j (2) 



( 3 ) 
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When smoke propagation in the space needs to be considered, an additional 
equation similar to (3) for the smoke concentration can be added (Yang and Liu, 
1978; Raycraft, 1987). It is, however, not shown here, since the smoke field is 
not considered in the demonstration of this field model, as will be shown 
later, and in the comparison of the numerical results and the test data. In 
Eqs. (1), (2) and (3), the physical coordinates are with i = 1, 2 and 3, 

p is the mean fluid density, u^. are the mean velocity components in x^, 
subscript t signifies the time derivatives, p is the mean static pressure, 

G^. represents the gravity vector components, T is the mean fluid temperature, 
and Q represents the source term. It is noted that inside the fire envelope, 

Q includes both the thermal radiation and combustion heat contributions, and 
outside the envelope, Q simply vanishes, provided that the gas radiation 
effects are neglected. Furthermore, the mean isobaric heat capacity Cp^, the 
shear stress tensor a.., and the conduction-flux q . are defined as 

1 J Cl 



c 



pm 




T 

c_(T)dT 

r 



( 4 ) 






( 5 ) 



XI 



= -k 



eff 



( 6 ) 



respectively, where is a convenient reference temperature generally taken 

to be the fluid temperature prior to the commencement of the fire, 

usual Kronecker delta, and and are the effective dynamic viscosity 
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and thermal conductivity, respectively, both accounting for the combined lami- 
nar and turbulent contributions. As will be shown later, both and 

need to be modeled separately. 

In addition, the fluid medium is taken to be a perfect gas, and the perfect 
gas relation 

P = pRT (7) 

where R is the fluid gas constant, is not only used to accommodate the 
compressibi.l i ty effects, but also needed to account for pressure build-up due to 
heating from the fire in a constant-volume space. As an approximation, all 
laminar transport properties, such as Cp^^, p and k of the fluid, are considered 
to be temperature independent and calculated at T^. This is justifiable because 
Cpm normally only a weak function of temperature, while p and k only 
contribute to a minor degree to the respective effective values due to the much 
larger turbulent effects expected in the fire phenomena. The initial conditions 
are simply that there is no motion of the fluid, i.e., u^. = 0, and T = Tp. 

The boundary conditions are the usual no-slip conditions at the enclosure 
interior wall surfaces and that for the temperature field is determined by a 
heat balance there. When the interior wall is heated up, heat loss occurs 
through the wall thickness and eventually to the environment from the exterior 
wall. In the present study, such wall losses are accommmodated by transient con- 
duction through the wall and then by combined convection and radiation at the 
outer surface by means of a constant coefficient of heat transfer. Since the 
wall thickness is always small compared to the characteristic length along the 
wall, a one-dimensional conduction analysis should be sufficient. Thus, for 
conduction through the wall thickness, we may write 



where the subscript refers to the wall material, and 
coordinate, together with the boundary conditions; 



n is the outward normal 



n = 0 



n 




8T 



9T. 



q - k = -k — : 

r eff 8n s an 



3T 



^ = h(T - Te) 



(9) 

( 10 ) 



where n = 0 is at the inner surface, L is the wall thickness, h is the 

w 

coefficient of heat transfer at the outer surface, is the environment tem- 
peratures and is the thermal radiation flux arriving at the surface from 
the enclosure interior. As pointed out previously, h is taken to be a con- 
stant, depending on the wind conditions outside. Before Eqs. (1), (2), (3) 
and (8) can be solved simultaneously, additional turbulence and thermal radia- 
tion models must be introduced, and these models will be given in a later sec- 
tion. 

Now the above governing equations in are ready to be transferred 

into those of a generalized orthogonal coordinate system given by the coor- 
dinates The two coordinate systems are related by means of the scale 

h^ in the direction 0^ given by 



factors 
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where is the base vector. One of the special properties of the orthogonal 
coordinates is that these base vectors are orthogonal so that the covariant and 
contravariant metric tensors satisfy the following conditions 

9ij - 9i • ij - 6ijhihj (12) 



gij 



6 . . 



(13) 



resulting in a diagonalized metric tensor, where g 
or 



is the determinant of 




g - igijl 






(14) 



With the above orthogonal properties , transformations of Eqs. (1), (2) and 
(3) in the system lead to the following respective equations in the 9^ 

system: 






(15) 



uV/h.) - - P_./h. + c<;i + ^ 



1 . . . . 3h. 



( 16 ) 



-9- 




p C T/h.) 
pm 1 




( 17 ) 



where 




(18) 



The above conservation equations are clearly valid for any orthogonal coor- 
dinate systems. In particular, h^ = h 2 = h^ = 1 for the Cartesian coor- 
dinates; h^ = r = 0*, h 2 = hg = 1 for the cylindrical coordinate system where 

r is the radial variable; and h^ = rsinej) = e^sina^, h^ = 1 , h^ = r = 9* for 

the spherical coordinates, where <p is the polar angle. The boundary and ini- 
tial conditions need no further transformation. In addition, in view of the 
relatively thin wall thickness, one-dimensional conduction is still valid and 
therefore, Eqs. (8), (9) and (10) remain the same. 

As pointed out previously, before the governing equations (15), (16), (17) 

and (8) can be solved, physical models for the radiation exchange in the enclo- 

sure and for the turbulence field are still needed. Since the radiation 
exchange depends on the specific geometry of the enclosed space, the radiation 
model will be described later after the geometry of the NRL FIRE I test facility 
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is introduced. However, a turbulence model utilized in the present study can b 
given here, again in terms of an arbitrary orthogonal coordinate system. This 
model is based on a mixing-length formulation originally developed for two- 
dimensional recirculating flows (Nee and Liu, 1978), and has been extended to 
three-dimensional flows with good success (Yang and Lloyd, 1985). In terms of 
the orthogonal coordinate system, the equation for the effective viscosity 
according to this model can be written as 



eff 



= 1 + 



n_. 1 n 



‘i 99- 



(19) 



2 + 



Ri 

Pr^ 



where p- is the laminar viscosity at T„, Pr is the turbulent Prandtl 



number taken to be unity, and R^. is the gradient Richardson number given by 



Ri = -4 



JI 









hn‘ 



3n 



( 20 ) 



where U„ is a reference veloci ty, and n is a unit vector in the negative 
gravity direction. In addition, the mixing length 1 normalized with the 
enclosure height H can be written as: 




(uV)V2 



1 3u^.2,l/2 



i,j j 36^ 

TT 



X,3 D 99^ i,j 1 j 36 36^ 



( 21 ) 



where K is a constant, normally taken to be 0.2. Finally, the effective 



eff 



conductivity k 



is given by 
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k 



eff 



1- + JL 



( 22 ) 



k 



Pr Pr^ Ur 



R 



where k„ is the molecular thermal conductivity at T_, and Pr is the 
R R 

corresponding laminar Prandtl number. 

NUMERICAL FINITE-DIFFERENCE CALCULATIONS 

In the present study, the governing equations (15), (16) and (17) in the 
generalized orthogonal coordinates are solved by means of the primitive-variable 
staggered control -vol ume finite-difference approach. As recently pointed out by 
Yang et al. (1987b, 1988), one major advantage of utilizing the fully transformed 
equations is that by judiciously choosing new definitions of flux and source 
quantities, the final finite-difference equations appear very much like those 
for the well-known Cartesian-coordinate cases, and hence can be handled in a 
routine manner for their solutions. That this is the case can be illustrated as 
fol lows; 

For a control volume given by AV = /g A0’A6®A0®, where A indicates a 
step size, the mom entum equation (16) for u’, after integrated over the 
staggered control volume surrounding a grid point P, becomes 




12 13 13 

- Ms'^As + M. A. - M. “^A. = S 
f f D b 



( 23 ) 



where the areas A are given by 
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e ,w 



f.b 



A 



n,s 



(hiA9^2Ae2)^_b 

(h3A0^h,A9^) 



( 24 ) 



n,s 



and is the total momentum flux along the 0 ^ direction for the velocity 

component u^ due to convection and diffusion. The subscripts e, w, n, s, f 
and b denote the east, west, north, south, front and back boundaries one-half 
cell size away from the point P, respectively. More specifically. 



= ( p - o^) 



(25) 



S = -PgAe + Pv,A„ + PG^AV - m12 (a„ - Ag) 



- Ml3(Af - Afa) + (m 22 + M33)(Ag - A^) 



(26) 



Here the idea of Raithby et al. (1986) of using a stress-flux formulation is 
utilized, and the momentum flux is then split into 



MID = + (o 3 -al) 

1 1 



(27) 



where 



a 



j 

i 



eff 



(h. (^)J 

3 90J 



M 



ij 



pu^u^ - 



3 

1 



(28) 



(29) 
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As a result, Eq. (23) now becomes 



where 



AV (PU) + + M^^A - M^^A 

P ee ww nH s£ 



^13 ^13 

* - .«b\ = s 



( 30 ) 



a = S - (o'l - - al)„A„ 






- (o, - ohj+ ( 






a, - a: 



'>b\ 



( 31 ) 



Here the stresses are evaluate(d from the information of the prior itera- 

tion, and the source term is known at the present iteration. Once these are 

properly evaluated, Eq. (30) becomes identical in form to that for the 
Cartesian coordinates. Further approximations to can be made. Basically, 

the convective terms are approximated by the QUICK scheme of Leonard (1983), 
extended to three-dimensional cases, and the diffusion terms are by central dif- 
ferences. The details have been given by Raycraft (1987), and hence will not be 
repeated here. 

The treatment of the energy equation is very similar. For example, a total 
heat flux j\ which includes both convection and conduction, is introduced as 



J 



i 




( 32 ) 



The corresponding finite-difference equation for the energy equation (17) then 



becomes 
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(p aV + Je Ab - Jw^ Aw + Jn^An 

- Js^As + J^^A^ - Jj^^Aj^ = Sj (33) 

where the source also includes Q and the solution procedure is the same 

as that for the momentum equations. 

In addition, it is noted that for a closed enclosure, any energy addition 
to the gas due to the fire tends to raise the pressure in the enclosure in view 
of the constant volume. Consequently, a time-dependent global pressure 
correction is needed to determine this rise, in addition to the local pressure 
correction in the primitive-variable control-volume procedure, which only 
depends on the local velocity variations. For this purpose, the technique based 
on the perfect gas relation as developed by Nicolette et al. (1985) is utilized 
in the present study, and this technique has also been described by Raycrdft 
(1987). 

NRL FIRE I TEST FACILITY AND RADIATION MODEL 

As pointed out previously, the validation of the field model developed in 
the present study is based on test data obtained in the Naval Research 
Laboratory (NRL) FIRE I full-scale test facility (Alexander et al., 1982). As 
shown schematically in Fig. 1, the facility has a cylindrical midsection with a 
radius of 9.6 feet (2.93m) and a length of 27.4 feet (8.352m) and hemispherical 
end caps. Thus, it is in the form of a pressure vessel having a total volume of 
11,639 cubic feet (324m^). The instrumentation consists of pressure trans- 
ducers, thermocouples, and radiometers. There is also capability to measure 
smoke obscuration levels, gas composition and humidity. The use of circulation 
fans can also be included in a test to determine the effect of ventilation. 



FRAME RAYS 
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Pressure transducers are located at the north and south ends of the test 
chamber, while the temperatures are monitored with thermocouple arrays located | 
inside the spherical end caps as shown in Fig. 1. These are chromel -alumel j 
thermocouples with diameters of 0.12 mm and have ceramic insulation enclosed in 
304 stainless steel jackets 1.0 mm in diameter. As will be shown later, these j 
temperature readings provide the data for comparison with the numerical calcula-| 
tions based on the field model . 

In the field of fire modeling, the importance of thermal radiation contri- 
bution to the energy transfer and turbulent buoyant flow processes involved is 
generally well recognized. However, accommodation of the radiation effects in 
fire models is always difficult and requires solutions to the radiative 
transfer equation for complex geometries and nongray-gas, nonhomogeneous, and 
nonisothermal behaviors (Yang, 1986). The presence of soot and smoke further 
complicates the situation, even though some simplification here is possible. 
Therefore, any reasonable simplification in including the radiation effects is 
desirable. In the present study, only surface-surface and flame-surface 
radiation exchange is includod in the field model. In other words, the gas out- 
side the flame envelope is taken to be transparent. In addition, the flame is 
considered to be represented by an equivalent gray surface which exchanges 
radiation with all other surfaces in the enclosure. This simplification is 
deemed reasonable, since it is known that the effects of participating gases in 
an enclosure is largely to equalize the gas temperatures, and the tdcliation 
exchange is in general dominated by surface-surface interactions (Chang et al., 
1983). Under this simplification, the radiation heat transfer can be readily 
calculated by any one of the well-known enclosure procedures (Sparrow and 
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Cess, 1978; Siegel and Howell, 1981). In the present study, the radiation 
calculations are based on the net radiosity method which is particularly suited 
for gray surfaces. The radiosity at a surface i can be written as 



B. 

1 



e 




+ 0 






N 

z 

j=l 



B.F. . 
J ij 



( 34 ) 



where e is the emissivity, and F. . is the direct view factor from surface i 
to surface j. The index N denotes the total number of surfaces in the enclo- 
sure, including that of the flame. The above equation results in N linear, 
inhomogeneous, algebraic equations for the N unknown radiosities. By solving 
the simultaneous equations, usually by matrix inversions, B^ can be found, 
and the radiation heat flux at surface i can then be readily calculated from 

e . 



- — (oT - B.) 



Vi 1 - e. ' i i 

1 



( 35 ) 



which, for all enclosure interior surface elements, becomes part of the thermal 
boundary conditions in accordance with Eq. (9). This method of calculating q^^ 
can be further simplified by obtaining a direct relation between q^^ and the 
temperatures T^ so that matrix manipulations can be kept at a minimum. Eq. 
(34) can be rewritten as 



where 




X. .B. 
iJ J 






(1 - e.)F. . 
i U 



S- 



( 36 ) 



( 37 ) 
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which, in a matrix form, becomes 



[X] <B> = o<T'*; 



( 38 ) 



By inverting the coefficient matrix [X], we obtain for the radiosity 



<B> = a [X]”l <t4> 



( 39 ) 



Then in accordance with Eq. (35), we finally have 



N 



^ G . ,qT . 
'■1 j=i ij j 



( 40 ) 



where 









-1 



( 41 ) 



Here ’ only depends on the emissivities which are taken to be constants. 



and consequently, together with can be calculated once for all in a given 

problem. 

It is thus clearly seen that the radiation model utilized in the present 
study essentially reduces to a problem of determining all the appropriate view 



factors F... To be consistent to the finite-difference calculations for the 
TJ 



buoyant flow in the enclosure, in which each interior surface cell has a tem- 
perature which depends on time, these surface cells can also be conveniently 
used as surfaces i for radiation. Consequently, the determination of the view 
factors does not require any surface integration, and they can be calculated 
directly from the definition 
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iA.Aj 



( 42 ) 



where /3 is the angle between the normal of the surface and R, the line con- 



enclosure, no generic forms can be written in terms of the generalized ortho- 
gonala coordinates, and this is the very reason that the radiation model cannot 
be introduced until the geometry is specified. 



cell pairs as internal cylindrical to cylindrical, cylindrical to hemispherical, 
and hemispherical to hemispherical surfaces. For the last case, surfaces i 
and j may also be located on different end caps. Surprisingly, not all such 
view factors are known in the literature (Siegel and Howell, 1981). As a 
result, these view factors have been derived directly from the FIRE I geometry. 
They have been given by Raycraft (1987) and there is no need to repeat them 
here. Furthermore, the view factors associated with the flame must be treated 
somewhat differently, and this will be given in the following section. 

RESULTS OF NUMERICAL CALCULATIONS AND DISCUSSIONS 

Detailed numerical calculations based on the field model have been carried 
out for the NRL FIRE I facility geometry so that the results can be compared 
directly with the test data. One set of test data for burning methanol inside 
the FIRE I tank is chosen for this purpose. For this burn test, the tank is 
essentially empty and the fuel pan is located at 23.1 feet (7.041m) from either 
end cap and elevated 3.21 feet (0.978m) from the bottom. The initial tem- 
perature and pressure are at 35.6®C and 1 atm., respectively. Shown in Fig. 2 



necting A. 



and A . . 
J 




For the geometry of FIRE I, the view factors F^j depend on such surface 
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Figure 2(a). Front View of Calculation Grid for FIRE I 




- 21 - 




Figure 2(b). End View of Calculation Grid for FIRE I 
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is the grid system utilized in the computations. As seen, there are 20 
segments in the circumferential direction, 12 radial segments representing the 
interior of the tank, and 18 segments along the cylindrical midsection, for a 
total of 6,720 interior cells. Also, in order to avoid calculation difficulties 
surrounding the line of zero radius, a row of cells is added to preserve 
the direction of the velocity vectors (Yang et al . , 1987b). The thickness of 
the wall is 0.375 in (0.95 cm), and occupies one single cell segment. 
Consequently, the entire interior surface of the tank has 560 surface cells 
which all exchange radiation with the flame and also with each other, and there 
are also 560 cells for conduction through the wall thickness. Furthermore, 
similar to that in the UNDSAFE-II field model (Yang and Liu, 1978), the flame is 
modeled as a volumetric heat source, and in the present calculations, the entire 
heat release rate in the flame is taken to be distributed uniformly over 19 
cells above the fuel pan. Thus, radiation heat transfer calcula-tions involve 
altogether 579 surfaces. The numerical calculations also require several physi- 
cal data relative to the wall properties and external boundary conditions. In 
addition to the wall thickness given previously, the specific heat, thermal con- 
ductivity, and density of the metallic wall (ASTM-285 Grade C steel) are taken 
to be 0.1 B/(1bF) (0.42 kJ/kgK), 25 B/hr.ft.F(43.27 W/mK), and 487 Ib/ft^ 

(7,801 kg/m’), respectively. There is also convection and radiation heat 
transfer occurring at the tank exterior surface, and the corresponding coef- 
ficient of heat transfer is taken to be a constant of 15 B/hr.ft®F (85.169 
W/m’K). This coefficient can be adjusted in the computations, if needed. The 
time increment is essentially determined by numerical stability requirement. 
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The driving force in the physical problem is the heat release rate from the 
fire. Since in the present model the flame is taken to be a volumetric heat 
source, the unsteady heat release rate must be treated as an input. Normally, 
such data can je obtained from the fuel depletion rate and the heating value of 
the fuel. Unfortunately, because of instrumentation error during the methanol 
fire test, the fuel depletion rate data were also in error, and consequently the 
heat release rate must be determined in an indirect way. Since the static 
pressure in the tank during the fire test was recorded, and in a constant volume 
heat addition process, the heat added is directly proportional to the rate of 
rising static pressure, an estimate of the heat release rate can be obtained 
from the recorded pressure data. In an attempt to generate the unsteady heat 
release rate in this way automatically on the computer, a pressure-tracking 
routine based on feedback control and the field model has been developed 
(Nies, 1986; Raycraft, 1987). The result is shown in Fig. 3. Even though the 
routine is performing well insofar as the pressure level is concerned, the 
resulting heat release rate, which is sensitive to the time derivative of the 
pressure, oscillates with large variations leading also to large variations in 
the temperature field in the tank. This heat release rate from the pressure- 
tracking is given in Fig. 4. Since such large variations of the heat release 
rate are not expected in the fire test, it is reasonable to use a curve-fitted 
variation as the input to the field model calculations, as also shown in Fig. 4. 
In order to demonstrate that such a curve-fitted heat release rate variation is 
a reasonable approximation to the real data, the calculated pressure rise is 
also given in Fig. 3, and it is seen that it does follow the experimental curve 
rather well. Detailed intermediate results based on pressure tracking can be 
found in the thesis of Raycraft (1987). 



PRESSURE ( ATM ) 
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Figure 3. Pressure History After Commencement of Fire 
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Figure 4. Input Heat Release Rate 



-26- 



A more important result of the numerical computations is to compare the 
temperatures at the various thermocouple locations with those from the test. 

The thermocouples are all mounted on vertical racks in the two end-cap regions 
close to the edges of the cylindrical midsection (Fig. 1). The three ther- 
mocouples, 1, 2 and 4, are located in the south hemispherical end cap, and the 
readings there are compared with the numerical results. Thermocouple 1 is 
located 79 inches (2.01 m) above the midplane of the tank. Thermocouple 2 is one 
foot (0.305 m) below thermocouple 1, and thermocouple 4 is at two feet (0.61 m) 
below thermocouple 2. The thermocouple 3 was defective during the test, and 
hence comparisons there are not made. Figs. 5, 6 and 7 show the temperature 
variations at thermocouples 1, 2 and 4 locations, respectively, as compared to 
those from the test, which represent smoothed mean temperatures. The agreement 
is reasonable at the two lower thermocouple locations with maximum deviations of 
less than 10*C. However, the calculations definitely overpredict the tem- 
peratures at the thermocouple 1 location by a maximum of about 30®C. Several 
possible reasons may be responsible. One is that since the thermocouple 1 loca- 
tion is closest to the tank wall, its temperature is most affected by the boun- 
dary condition at the tank exterior surface in terms of outside ambient 
temperature and the coefficient of heat transfer, both of which have been estimatl 
calculation purposes. Another possible reason is that in the model, the heat 
release rate is taken to be uniformly distributed. In reality, considerable 
vertical nonuniformity may result from the combustion process in the flame. If 
more heat is released in the lower pari: of the flame, more heat would be 
radiated away there so that the thermal plume would attain a lower temperature, 
resulting also in lower temperatures at the thermocouple 1 location. A third 
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Figure 5. Temperature Comparison at Thermocouple #1 Location 
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Figure 6. Temperature Comparison at Thermocouple #2 Location. 



120 



-29- 




O 

(D 



O 



O 

C\J 



O 

O 



O 

00 



o 

(£> 



O 



O 

C\J 



O 



CO 



U 

c 

o 

o 

CD 

CO 




( do ) 3afUVU3dW3i 



Figure 7. Temperature Comparison at Thermocouple #4 Location 
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possible reason is that gas and soot radiation is neglected in the field model. 
Since in general the effect of gas and soot radiation is to make the tempera- 
ture field more uniform in the enclosure (Yang, 1986), the thermocouple 1 loca- 
tion would also tend to have lovj* ; lemper atures. Also, the uncertainty in the 
experimental heat release rate itself could also be responsible for the degree 
of discrepancy shown in Fig. 5. Incidentally, it is noted in Figs. 5, 6 and 7 
that all calculated temperatures exhibit slow oscillations which cannot be 
attributed to turbulence. Such oscillations could very well be real because 
of inherent thermal instability .associated with thermal plumes (Yang et al., 
1983) . 

One of the major advantages of the field model is the fact that the numeri- 
cal solution gives detailed spatial and temperal variations of the velocity 
and temperature fields. Such information is critical in the development of 
fire mitigation measures, especially in the placement of barriers and par- 
titions, and in the determination of fire loads on the walls. The field 
results from the calculations in terms of the velocities and temperaturess can 
now be shown. Because of the three-dimensional nature of the fields, it is 
necessary to look at these quantities at specific sections of the tank. At a 
time instant 30 seconds into the fire, the isotherms, or the temperature field, 
and the velocity vector field are shown in Figs. 8. Fig. 8(a) gives a front 
view of these fields at the mid section of the tank. Fig. 8(b) refers to the 
view also at the mid section from the south end of the tank, while Fig. 8(c) 
is another view from the south end, but at the intersection between the end 
cap and the cylindrical mid section, or at the base of the end cap. Some 
sections are used in Figs. 9, 10, 11 and 12 at 60 seconds, 90 seconds, 120 
seconds and 150 seconds from the coinmoncement of the fire, respectively. 
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Figure 8(a) . 






Mid-Section Front View of Isotherms and Velocity 
Field at 30 seconds 




! 





Figure 8(b). Mid-Section End View of Isotherms and Velocity Field at 30 seconds 
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Figure 8(c). Section View at Base of End Cap 
at 30 seconds 




Figure 9(a). flid-Section Front View of Isotherms and Velocity 
Field at 60 seconds. 
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Figure 9(b). Mid-Section End View of Isotherms and Velocity Field at 60 seconds 
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Figure 9(c). Section View at Base of End Cap of Isotherms and Velocity Field at 60 seconds 
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Figure 10(a). Mid-Section Front View of Isotherms and Velocity 
Field at 90 seconds 
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Figure 10(b). Mid-Section End View of Isotherms and Velocity Field at 90 seconds 
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Figure 11(a). 



Mid-Section Front View of Isotherms and Velocity Field 
at 120 seconds. 
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Figure 11(b). Mid-Section End View of Isotherms and Velocity Field at 
120 seconds. 
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Figure 11(c), Section View at Base of End Cap of Isotherms and Velocity Field at 120 seconds 
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Figure 1 2(a) . 



Mid-Section Front View of Isotherms and Velocity 
Field at 150 seconds 
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Figure 12(b). Mid-Section End View of Isotherms and Velocity Field at 150 seconds 
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150 seconds. 
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Since the fire is located at the center of the tank, both sections in 
Figs. 8(a) and 8(b) include the flame, and therefore it is seen that the fire 
plume is already well formed at 30 seconds. There is also a hot ceiling layer 
along the ceiling of the tank, and the lower two thirds of the tank are still 
at almost the initial temperature. The general flow circulation pattern is 
much more complex than might be expected. The entrainment into the fire plume 
is limited to the immediate neighborhood of the fire, as can be seen in Figs. 
8(a) and 8(b). The flow in the hot ceiling layer does not seem to have strong 
enough momentum to carry the flow toward the lower half of the tank, but re- 
turns into the tank interior, resulting in a downward-biased flow even at the 
section as shown in Fig. 8(c). The flow completes its journey back toward the 
fire region in a somewhat spatially oscillatory path. At 60 seconds as shown 
in Fig. 9, the depth of the hot ceiling layer is increased, even though the 
lower half of the tank is still not affected much by the fire. The flow in 
the ceiling layer toward and down the spherical end cap is definitely spiral 
in nature. In addition, the spatially oscillatory flow back to the fire zone 
become sufficiently large that it induces a backward flow toward the end cap 
again at the bottom of the fire. It may also be of interest to note that the 
isotherms are now more densely packed at the ceiling region of the tank, indi- 
cating that heat losses through the wall now become more important. Similar 
trends of the temperature and velocity fields can be seen at 90, 120 and 150 
second instants in Figs. 10, 11 and 12, respectively. It is interesting to 
note that the hot ceiling layer reaches down almost to the half tank height 
level at 120 and 150 seconds, and inside the hot layer, the temperature field 
is essentially stably stratified except a small thermal bounday layer region 



-47- 



next to the ceiling. Incidentally, the isotherm and velocity vector scales 
in Figs. 8 to 12 are not uniform, and consequently the diagrams should not be 
used quantitatively for comparison purposes. 

CONCLUDING REMARKS 

The present study develops a three-dimensional differential field model based 
on general orthogonal coordinate systems which includes the physical effects of 
turbulence, strong buoyancy, full compressibility, surface-surface and surface- 
flame radiation exchange, pressure rise due to fire loading, and heat losses 
through the wall. The model is capable of dealing with geoometries that 
represent enclosures with combinations of shapes which are describable by 
orthogonal coordinates. The model, which is also based on a control -volume 
finite-difference approach with primitive variables, staggered cells, and QUICK 
scheme, and both local and global pressure corrections, is used to simulate a 
set of full-scale fire test data obtained in the NRL FIRE I fire test facility, 
which is a pressure vessel with a cylindrical mid-section and hemispherical end 
caps on both ends. Despite some uncertainty in the heat release-rate input 
data, reasonable comparisons in both the temperature levels at selected points 
in the test tank and the tank pressure between the calculated and test results 
have been obtained. However, there is also indication that the calculated tem- 
peratures close to the tank ceiling overpredict the test values, and possible 
reasons for this discrepancy have been noted. 

Since the completion of this study, two additional capabilities of the 
field model have been developed and these include effects of recirculating 
fans to provide ventilation inside the test tank and the placement of 
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horizontal deck inside the tank. At the present time, the field model treats 
the fire as a volumetric heat source with prescribed heat release rate. 
Efforts have been initiated to develop a combustion model so that the fire 
envelope and heat release rate do not need to be prescribed. Furthermore, 
plans are also being made to incorporate a gas and soot (smoke) radiation 
model and a more refined turbulence model into the overall field model. 
Results of these efforts will be reported at future dates. 
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NOMENCLATURE 



A 


boundary areas of control volumes, m® 


A. 

1 


radiating surface, m^ 


B. 


radiosity, w/m® 


c 

P 


isobaric heat capacity, kJ/kgK 


c 

pm 


mean isobaric heat capacity, Eq. (4), kJ/kgK 


F. . 
ij 


radiation view factor 


G . 
1 


components of gravity vector, m/s^ 


G. . 


radiation matrix, Eq. (41) 


g 


determinant of g . . 

ij 




covariant metric tensor 
contravariant metric tensor 


g 


gravity vector 




base vector 


H 


enclosure height, m 


h 


coefficient of heat transfer, w/m*K 


h . 
1 


scale factor 




total heat flux, w/m^ 


K 


constant in turbulence model 


k 


thermal conductivity, W/mK 


L 

w 


thickness of wall, m 


! 


mixing length, Eq. (21), m 


M^J 


total momentum flux, N/m^ 


M^J 


momentum flux, Eq. (29), N/m^ 


n 


outward normal 


Pr 


Prandtl number 
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P 

Q 



R . 

1 

R . . 
T J 

r 

S 



u . 

1 



X . . 

x:^. 

ij 



X . 
1 



|3 



1 



6 , . 

T J 



€ . 
1 



0 

1 

0 ' 
p 



turbulent Prandtl number 

static pressure, N/m^ 

heat source, w/m^ 

conduction flux, w/m^ 

radiation flux, w/m’ 

gas constant, Nm/kgK 

gradient Richardson number, Eq. (20) 

line between surfaces i and j 

radius variable, m 

source term, N 

source term, Eq. (31),N 

absolute temperature, K 

outside ambient temperature, K 

velocity components in X^. system, m/s 

velocity components in 0^ system, m/s 

matrix components, Eq. (37) 

i nverse of X . . 

Cartesian coordinates, m 

angle between normal to A. and R.., rad 



1 J 



control volume, m^ 

Kronecke delta 
surface emissivity 
circumferential angle, rad 
orthogonal coordinate 
viscosity, kg/ms 
density, kg/m’ 
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a 

a. . 



ij 



a^. 



ij 



aJ 
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Subscripts 

eff 



R 



s 

t 



Stef an-Bol tzmann constant, w/m^K'* 
shear stress tensor, N/m^ 
shear stress tensor ! r 0 system, N/m® 
shear stress tensor, Eq. (28), N/m’ 
polar angle, rad 

effective (laminar + turbulent) 

reference quantities 

solid 

time derivative 
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